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In  this  report  we  investigate  the  accuracy  with  which  the 
position  of  a  receiver  can  be  determined  by  use  of  a  high-altitude 
satellite  navigation  system.  The  navigation  system  is  modelled 
as  a  problem  of  nonlinear  estimation  in  the  presence  of  random 
disturbances.  Equations  are  derived  for  describing  positioning 
errors  by  using  linearization  and  the  Kalman-Bucy  filtering 
equations. 


Accepted  for  the  Air  Force 

Franklin  C.  Hudson 

Chief,  Lincoln  Laboratory  Office 


Table  of  Contents 


1 .  Introduction  1 

Z.  Linear  Estimation  Theory  3 

Z.  1  The  Linear  Estimation  Model  3 

Z.  Z  The  Linear  Estimation  Problem  4 

Z.  Z.  1  Special  Case  I:  Time- Independent  Observations  6 

Z.  Z.  Z  Special  Case  II:  Time -Dependent  Observations  8 

Z.  Z.  3  General  Case:  Time-Dependent  Observations 

and  Parameters  1  1 

3.  Nonlinear  Estimation  Theory  16 

3.  1  The  Nonlinear  Estimation  Model  16 

3.  Z  The  Nonlinear  Estimation  Problem  16 

3.  Z.  1  Special  Case:  Time-Independent  Observations  19 

3.  Z.  Z  General  Case:  Time -Dependent  Observations  21 

4.  Application  of  Nonlinear  Estimation  Theory  to  the  Problem  of 

Navigating  with  High- Altitude  Satellites  Z4 

4.  1  Special  Case  I:  Fixed  - Po sition  and  Single  Observation  24 

4.  Z  Special  Case  II:  Simple  Motion  and  Multiple  Observations  Z9 

Appendix  33 

References  34 


v 


NAVIGATION  WITH  HIGH- ALTITUDE  SATELLITES: 

A  STUDY  OF  ERRORS  IN  POSITION  DETERMINATION 

1.  INTRODUCTION 

A  possible  scheme  using  high -altitude  satellites  for  very  accurate 
position  determination  and  navigation  on  the  earth’s  surface  has  been  des¬ 
cribed  by  Goblick.  ^  With  this  scheme,  several  satellites  independently  trans¬ 
mit  timing  signals  generated  from  free-running  clocks.  Control  stations 
monitor  the  transmission  of  timing  signals  from  the  satellites  and  compare 
the  observed  timing  signals  to  an  accurate  master  clock.  Satellite  clock 
timing  errors  are  then  transmitted  to  users  by  the  control  stations  via  the 
navigation  satellites  themselves  and  in  the  form  of  digital  data  superimposed 
on  the  navigation  signals.  The  satellite  clocks  are  thereby  effectively  syn¬ 
chronized  by  the  control  stations.  Moreover,  the  control  stations  monitor 
the  satellite  positions  and  also  relay  this  information  to  users  as  digital  data. 
A  user  with  a  synchronized  clock  could  determine  his  position  by  knowing*. 

(a)  his  distance  from  the  earth's  center;  (b)  the  satellite  positions;  and  (c)  the 
arrival  time  of  the  timing  signals  from  two  satellites.  If  the  user  does  not 
have  a  synchronized  clock,  three  satellites  are  required  and  differences 
between  the  arrival  times  of  timing  signals  from  two  pairs  of  satellites  are 
used . 

The  scheme  has  the  advantage  of  requiring  only  passive  operations  by 
the  user.  In  addition,  position  fixes  can  be  accomplished  in  a  relatively 
short  period  of  time  compared  to  systems  employing  a  low-altitude  satellite. 
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In  this  report  we  shall  investigate  the  accuracy  with  which  a  user's 

position  can  be  determined  in  the  presence  of  uncertainties  which  inevitably 

exist;  the  uncertainties  arise  both  because  of  observation  noise  which,  for 

example,  results  in  imprecise  time  or  range  measurements  and  because  of 

inaccurate  knowledge  of  the  satellite  positions,  user's  height,  and  timing. 

We  consider  two  cases.  In  the  first,  we  assume  the  satellites  and  user  are 

motionless  and  that  only  one  observation  is  made.  *  For  the  second  case,  we 

account  for  simple  satellite  and  user  motion  as  well  as  multiple  observations. 

In  addition  to  studying  positioning  errors,  we  also  indicate  the  structure  of 

processors  for  realizing  the  position  fix. 

The  theory  we  shall  present  incorporates  a  linearized  version  of  the 

navigation  problem  so  that  it  is  convenient  to  divide  the  discussion  into  three 

steps.  In  the  first  step,  we  consider  the  general  linear  estimation  problem 

formulated  in  a  way  that  results  in  the  discrete  filtering  equations  of  Kalman 
2-3 

and  Bucy  .  The  derivation  of  these  equations  has  been  given  by  several 
2 

authors  ,  so  for  brevity  we  present  only  the  broad  points  leaving  certain 
of  the  mathematical  subtleties  to  the  cited  references.  In  the  second  step, 
we  consider  the  general  nonlinear  estimation  problem  by  a  technique  employ¬ 
ing  linearization  and  the  application  of  Kalman  and  Bucy's  equations.  This 

7 

technique  appears  to  have  first  been  used  by  Smith  et  al.  More  recent 

9  10 

discussions  are  given  by  Ohap  and  Cox  .  Finally,  in  the  third  step,  we 
show  how  the  navigation  problem  fits  into  the  model  of  the  general  nonlinear 
estimation  problem  and  then  simply  apply  the  previously  derived  results. 

*  This  portion  of  the  study  is  based  on  unpublished  studies  performed  at  MIT 
Lincoln  Laboratory.  See  the  acknowledgment  on  page  35. 
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2.  LINEAR  ESTIMATION  THEORY 


We  now  postulate  a  model  for  a  general  linear  estimation  problem.  The 
interpretation  of  the  model  is  left  to  Section  4  where  the  navigation  problem 
is  presented.  For  the  present,  we  simply  mention  that  the  model  is  a  state- 
variable  representation  which  can  describe  a  wide  class  of  discrete  linear 
systems  and  their  observed  responses  to  Gaussian  input  sequences. 

2.  1  The  Linear  Estimation  Model 

Let  x(i)  and  r(i)  for  i  =  1,  2,  .  .  .  ,  m  be  vector  sequences  of  Gaussian 
random  variables  defined  by: 

x(i)  =  0(i,  i-l)x(i-l)  +  G(i)  w(i)  (1) 

x(0)  =  xQ 

and 

r(i)  =  H(i)x(i)  +  n(i)  (2) 


where,  for  i  =  1,2,...  ,  m: 

(a)  x(i)  is  an  k-dimensional  column  vector.  The  initial  value  for  the 

sequence  is  a  random  variable  denoted  by  x^.  We  assume  that  the  a  priori 
knowledge  about  consists  of  a  mean,  E[x^]  *  x(0),  and  covariance  matrix, 
E[x0-x(°)]  [x0-x(0)]  4  V(0); 

(b)  ^(i,  i-  1)  is  a  known  transition  matrix  relating  x(i)  to  x(i-  1): 

(c)  w(i)  is  an  X -dimensional  uncorrelated  vector  - Gaus  sian  sequence 
with  zero  mean  and  known  covariance  matrix: 


E[w(i)w,(j)]  =  W(i)6..  (3) 

where  W(i)  is  a  symmetric  non-negative  definite  IX  X  matrix  and  is 
unity  when  i=j  and  zero  otherwise.  Deterministic  and  correlated  input 
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sequences  can  be  accommodated  by  making  w(i)  have  a  non-zero  mean  and 
by  increasing  the  dimension  of  x(i),  respectively; 

(d)  G_(i)  is  a  known  k  X  1  matrix; 

(e)  r(i)  is  an  observed  p-dimensional  vector  Gaussian  sequence.  We 
shall  denote  the  collection  of  observed  vectors,  r(l),  r(2),  .  .  .  ,  r(m),  by 

-1,  m; 

(f)  H(i)  is  a  known  p  Xk  matrix; 

(g)  n(i)  is  a  p-dimensional  vector  Gaussian  sequence  with  zero  mean 
and  known  covariance  matrix: 

E[n(i)rf(j)]  =  N(i)  6^  (4) 

where  N(i)  is  a  symmetric  positive  definite  p  X  p  matrix;  n(i)  represents 
noise  interferring  with  the  observation  of  H(i)x(i).  By  assuming  N(i)  to  be 
positive  definite,  we  are  implying  that  no  disturbance -free  observations  can 
be  made. 

2.  2  The  Linear  Estimation  Problem 

The  problem  we  wish  to  consider  is  that  of  estimating  x(m),  the  state 

th 

of  the  linear  system  at  the  m  point  in  the  sequence,  *  given  a  specified 
criterion  which  the  estimate  is  to  satisfy  and  the  m  observations 
[r(i):  1  ^  i  ^  m]  =  r^  We  shall  use  a  maximum  a  posteriori  probability 

criterion- -the  estimate  is  that  value  of  x  maximizing  p(x  |  r^  However, 

because  of  the  linearity  and  Gaussian  assumptions,  our  estimate  is  identical 
to  that  resulting  from  other  criteria  such  as  minimum  mean  square  error. 

*  As  new  data  is  observed,  m  increases  in  real  time. 
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Before  proceeding,  it  is  convenient  to  give  some  indication  of  the 
approach  to  be  taken  by  first  considering  the  following  two  special  cases  of 
the  general  problem  of  interest. 

Case  I:  For  the  first  problem,  we  remove  all  time  dependence  and  let 

r  =  Hx  +  n  (5) 

where  x  and  n  are  uncorrelated  zero-mean  Gaussian  random  variables  with 
positive  definite  covariance  matrices  W  and  N  respectively.  The  probability 
densities  of  x  and  n  are  then  given  by 

p(x)  =  - - - exp  [  —  t  x*  W.  1  x]  (6) 

V(2 n)m  det(W) 

and 

p(n)  =  —  1  -  exp  [  -  |n'  N  'n]  (7) 

^'(2tt)P  det(N) 


Case  II:  For  the  second  problem,  we  keep  x  fixed  and  reintroduce 
multiple  time-dependent  observations  by  letting 


r(i)  =  H(i)x  +  n(i)  i  =  1,  2,  .  .  .  ,  m 


(8) 


where  x  and  n(i),  i  =  1,2,...,  m,  are  uncorrelated  zero-mean  Gaussian 

fandom  variables  with  positive  definite  covariance  matrices  W  and  N(i), 

i  =  1,2,...,  m.  The  joint  probability  density  of  the  vectors  {n(i)  :  1,  2,  .  .  .  ,  mj 

n.  is  given  by: 

— l,  m  8  y 


p(n.  ) 
r  — l,  m 


. n  /( 2tt)P  det  [  N(i)] 
1=1  >  - 


exp 


-  i  Y  n'(i)N_1(i)n(i) 


i  =  l 


(9) 
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The  first  problem  is  illustrative  of  the  approach  and  is  also  of  interest 
because  of  its  close  relation  to  Schweppe's  error  analysis.  The  results 
from  this  problem  will  eventually  be  applied  to  the  navigation  problem  when 
both  the  satellites  and  user  are  fixed  and  only  one  observation  is  used  to 
determine  the  user's  position.  The  second  problem  is  presented  so  that  the 
notions  and  advantages  of  recursive  estimation  procedures  can  be  seen. 

2.  2.  1  Special  Case  I:  Time-Independent  Observations 

We  take  as  the  optimum  estimate  that  value  of  x  maximizing  the 
a  posteriori  probability,  p(x|r),  or  equivalently  its  logarithm,  In  p(x|r). 

Using  Baye's  rule,  we  have 

In  p(x  |  r)  =  In  p(r  |  x)  +  In  p(x)  +  cq 

where  c^  is  a  normalization  constant  independent  of  x.  From  Eqs.  6  and  7, 
we  then  have 

lnp(x|r)  =  —  y  [  r  —  Hx]  1  N  1  [  r  —  Hx]  —  yx'W  ^  x  +  c^'  (10) 

We  now  expand  the  right  side  of  Eq.  10  and  then  complete  the  square,  the 
result  being 

In  p(x|r)  =  “  \  [“  x'  H'  N' 1  r  +  x'  {  W "  1  +  H'  N"  1  H}  x  -  r*  N"  1  Hx]  +  c  " 
=  -  A  [x  -  VH'N'1  r]  1  V"1  [x  -  VH1  N~  1  r]  +  c"  (1/1) 

where  we  have  set 

V  =  {w_1  +  H' N'1  H}'1  (12) 

assuming  the  inverse  exists. 
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Inspection  of  Eq.  1  1  reveals  that  the  a  posteriori  density  of  x  is  normal 
and  maximized  by  the  conditional  mean,  VH'N  ^  r.  Consequently,  the  estimate 
we  seek  is  given  by 

X  =  VH'N_1r  (13) 

Furthermore,  the  error  covariance  matrix,  defined  as  the  conditional  expecta¬ 
tion  of  [x  —  x] [x  —  x]  ',  is  just  the  conditional  covariance,  V. 

It  is  observed  from  Eq.  12  that  the  evaluation  of  V  requires  the  inversion 
of  three  matrices.  An  alternative  expression  requiring  only  one  inversion 
can  be  given  by  applying  the  matrix  relation  proved  in  Appendix  1.  The  result 
is 

V  =  W  -  W  H'  [HWH1  +  N]  "  1  H  W  (14) 

Since  V  does  not  depend  on  r,  it  can  be  determined  before  making  an 
observation.  Consequently,  the  error  performance  can  be  studied  without 
actually  carrying  out  experiments.  Furthermore,  we  see  from  Eq.  13  that, 
if  V  is  determined  beforehand,  the  computation  of  x  involves  only  a  matrix 
multiplication  of  the  observed  data. 

Collecting  all  results,  we  have 

Model:  r  =  Hx  +  n 

Estimation  Equations: 

x  =  VH*  N_1r 

V  =  [W"1  +  H' N"1  H]  _1 

=  W  -  W  H'  [HWH'  +  N]  "  1  HW 
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Let 


e 


x  — 


A 

X 


be  the  error  vector  associated  with  x.  Then  the  mean- square  values  of  the 

components  of  e  are  the  diagonal  elements  of  V;  that  is,  Efe.^1  =  v...  The 
—  —  1  11 

error  vector  e  is  a  zero-mean  vector  Gaussian  random  variable  with  multi¬ 
dimensional  ellipsoids  as  constant  probability  contour s - -the  "error"  ellip¬ 
soids  are  defined  by 

e'  V  e  =  constant. 

The  squared  lengths  of  the  semiaxes  of  the  particular  ellipsoid, 
e'Ve  =  1, 


are  equal  to  the  eigenvalues  of  V  and  provide  a  convenient  coordinate -free 
measure  of  the  accuracy  in  estimating  x.  In  some  instances,  the  overall 
shape  and  orientation  of  the  error  ellipsoid  is  of  interest  and  can  be  examined 
by  determining  all  the  eigenvalues  and  eigenvectors  of  V,  respectively.  On 
the  other  hand,  simpler  measures  of  the  error  are  often  sufficient.  One  such 
measure  is  the  squared  length  of  the  semimajor  axis  which  equals  the  maxi¬ 
mum  eigenvalue  of  V;  this  measure  upper  bounds  Max  [v..]  .  Another  simple 

—  i  n 

measure  is  the  sum  of  the  squared  lengths  of  the  semiaxes  and  equals  the  trace 
m 


of  V, 


v. . . 

li 


Z.  Z.  Z  Special  Case  II:  Time-Dependent  Observations 

We  again  take  as  optimum  that  value  of  x  maximizing  the  a  posteriori 
probability,  p(x|r^  m),  or  equivalently  its  logarithm,  lnp(x|r^  ^).  Using 
Baye’s  rule  and  Eqs.  6  and  9,  we  have 
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lnp(x  jr^  m)  =  -  \ 


ill  <  I-. 

Yj  [r(i)  -  H(i)x]  1  N'1  (i)[r(i)  -  H(i)x]  -  £  x'  W-  x  +  cq 

L  —  1 

(15) 


where  Cq  is  a  normalization  constant  independent  of  x.  Proceeding  once 
again  through  the  steps  leading  from  Eq.  10  to  Eqs.  11  and  1Z,  we  obtain 
from  Eq.  1 5 


wh  ere 


(16) 


(17) 


It  can  be  easily  shown  that  x  is  an  unbiased  estimate  of  x.  V(m)  is  the 
associated  error  covariance  matrix;  that  is,  V(m)  is  the  conditional  expecta¬ 
tion  of  [x  —  x ( m ) ]  [x  —  x ( m ) ]  f. 

Equations  16  and  17  are  the  estimation  equations.  Note  that  all  the 
observations  {r(i):  1  <  i  <  m}  are  required  to  compute  x(m).  Data-storage 
requirements,  therefore,  increase  in  direct  proportion  to  the  number  of  observa¬ 
tions  taken.  Recursive  equations  for  £(m)  can  be  obtained  from  Eqs.  16  and 
17  by  simple  matrix  manipulation.  These  equations  provide  an  updating 
procedure  by  which  only  the  most  current  observation  is  used  in  conjunction 
with  the  preceding  estimate,  x(m-l),  and  covariance  matrix,  V(m-l),  to 
determine  x(m)  and  V(m).  They  have  the  advantage  that  storage  require¬ 
ments  do  not  increase  with  each  new  observation.  We  shall  derive  the  recursive 
equation  for  V(m)  first.  From  Eq.  17, 
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V(m)  =  [W'1  +  "V  _H'(i)  N_1(i)H(i)  +  H'(m)N_1(m)  H(m)l  _1 

1  =  1 

=  [V'1(m-l)  +  H,(m)N'1(m)H(m)]'1  (18) 

Equation  18  is  the  desired  relation.  An  alternative  form  can  be  obtained  by- 
use  of  the  lemma  of  Appendix  1,  the  result  being 

V(m)  =  V(m-l)  -  V(m-l)  H'(m)  [  H(m)  V(m- 1 )  H' (m)  +N(m)]  "  1  H(m)  V(m-l) 

(19) 

The  recursive  equation  for  x(m)  is  derived  starting  with  Eq.  16  from 
which  we  obtain 

x(m)  =  V(m)  H'(i)  N' ](i)  r(i) +  H'(m)  N' ](m)  r(m) 

Li=1 

=  V(m)  V"  1  (m  -  1 )  x  (m  -  1 )  +  V(m)  H'(m)  N~  1  (m)  r(m)  (20) 

From  Eq.  1 8  it  can  be  verified  that 

V(m)  V”  1  (m-  1 )  =  I  -  V(m)H'(m)N'1(m)H(m;  (21) 

Therefore, 

x(m)  =  x(m4)  +  V(in)H;(ffi)N'1(m)[i(m)-  HJm)x(m-l)]  (22) 

which  is  the  desired  result. 

Collecting  all  results,  we  have 

Model:  r(i)  =  H(i)x  +  n(i)  i  =  i,  Z,  .  .  .  ,  m 

Estimation  Equations: 

x(m)=x(m-l)  +  V(m)H,(m)N  ^  (m)  [  r(m)  —  H(rr. )  x(m-  1 )] 

V (m)  =  [V  ^(m- 1)  +  H'(m)  N  ^(m)H(m)]  ^ 

=  V(m-l)-  V(m-l)  H'(m)[H(m)  V(m-l)  H'(m)  +  N(m)]  _1  H(m)V(m-  1) 
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2.  2.  3  General  Case:  Time-Dependent  Observations  and  Parameters 

We  now  return  to  the  original  problem  and  derive  the  discrete  Kalman 
and  Bucy  equations.  In  the  following  procedure,  we  shall  derive  the  recursive 
equations  directly.  We  begin  by  observing  that 


p[  r(m),  x(m),  r 


— 1 ,  m  -  1 


p[r(m),  rlfin-1] 


pf  r(m)  lx(m),  r^  m_  j]  pfx(m)  |  m_1l 


p[r(m)!rUmll 

The  first  factor  in  the  numerator  can  be  simplified  when  it  is  noted  that  r(m) 
is  only  a  function  of  n(m)  when  x(m)  is  specified.  Then,  because  of  the 
independence  of  n(m),  x(i),  and  n(i)  for  i  <  m,  it  follows  that  r(m)  is  inde¬ 
pendent  of  r^  ^  i  when  x(m)  is  specified.  Consequently, 

p[r(m)  |x(m),  m_1]  =  p[r(m)|x(m)] 

and 


p[r(m)  lx(m)]  p[x(m)  |  r,  ,] 

p[x(m)jr.  ]  =  - = - - - , - = - ■  ’  m~L 

-1»  m  pfr(m)  |r^  m_1] 

We,  therefore,  have  the  result 


(23) 


In  p  [  x(m)  |  r.  ]  =  In  p  [  r  (m)  I  x(m)]  +lnp[x(m)|r.  ]  +c  (24) 

where  c^  is  a  normalization  constant.  The  first  term  on  the  right  is  given  by: 
In  p  [  r(m)  |x(m)]  =  —  ^  [  r  (m)  —  H(m)  x(m)]  1  N  ^(m)[r(  m)  —  H(m)x(m)]  +  c  1 


(25) 
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where  again  cj  is  a  normalization  constant.  We  now  need  to  investigate  the 

second  term  on  the  right  of  Eq.  24.  We  observe  first  that  x(m)  conditioned 

on  r.  .  is  normal  with  mean 
—1,  m-  1 

E[x(m)|r1  m_j]  =  x(m  |  m-  1 )  (26a) 

=  E[^(m,  m- 1)  x(m- 1)  +  G(m)^<m)  1^  m_j] 

=  0(m,  m- 1)  x(m- 1  |  m- 1)  (26b) 

where  x(m|m-l)  and  x(m-l|m-l)  denote  the  maximum  a  posteriori  estimates 
of  x(m)  and  x(m-l),  respectively,  based  on  the  m-1  observations,  r^  ^ 

In  deriving  Eq.  26,  we  have  used  the  independence  of  w(m)  and  rj  m  j  and 
the  identity  of  the  conditional  mean  and  maximum  a  posteriori  estimate.  The 
conditional  covariance  matrix  associated  with  x(m)  is  given  by 

F(m)  =  E[x(m)  —  x(m  |m- 1)]  [x(m)  —  x(m  |  m- 1 )]  1 

=  E  {^(m,  m-l)[x(m-l)  —  x(m  -  1  |  m  -  1 )]  +  G(m)  w(m) } 

X  {/(m,  m-l)[x(m-l)-  x(m-l  |m-l)]  +G(m)w(m)}’ 

=  0(m,  rn-l)V(m-l)^t  (m,  m-l)  +  G(m)W(m)G’(m)  (27) 

where  V(m-l)  is  the  error  covariance  matrix  defined  by 

V(j)  =  E[x(j)  -  x(j  |  j)]  [x(j)  -  x(j  |  j)]  '  (28) 

for  j  =  m  -  1 . 

Using  these  results  (Eqs.  25,  26,  and  27),  we  can  express  lnp[x(m)|rj  m]  * 
defined  by  Eq.  24,  as 
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lnp[x(m)|r.  ]  =  —  y[r(m)—  H(m)x(m)]'N  *  (m)  [  r(m)  —  H(m)  x(m)] 

—  y  fx(m)  —  x(m|m-l)]'  P  ^  ( m )  [  x  ( m )  —  x  ( m  |  m  - 1 )  ]  +  c^" 

(29) 

To  obtain  an  expression  for  x(m|m),  we  expand  Eq.  29  keeping  only  terms 
which  depend  on  x(m)  and  then  we  complete  the  square.  The  details  are 
straightforward  and  result  in  the  following  recursive  expressions  for  x(m|m) 
and  V(m): 

x(m  |  m)  =  V(m)[H'(m)  N_1  (m)  r(m)  +  P'1  (m)x(m|m-l)j  (30) 

and 

V(m)  =  [H’(m)  N"1  (m)  H(m)  +  P'^m)]"1  (31) 

where  P(m)  is  defined  by  Eq.  27. 

By  simple  matrix  manipulations,  Eqs.  30  and  31  can  be  placed  in  a 
form  having  some  computational  advantages  which  are  mentioned  below. 

Observe  from  Eq.  30  that 

x(m|m)  =  V(m)[Hf(m)N  ^(m)  r(m)  +  {H’(m)N  ^(m)H(m)  +  P  1  }  x  ( m  |  m  - 1 ) 

—  H’(m)N  \m)H(m)x(m|m-l)]  (32) 

=  x(m|m-l)+V(m)H,(m)N'1(m)[r(m)-  H(m)  x(m  |  m- 1 )]  (33) 

Equation  33  is  the  desired  expression  for  x(m|m).  By  using  the  lemma  of 
Appendix  1,  we  obtain  the  desired  expression  for  V(m)  from  Eq.  31 

V(m)  =  P(m)  —  P(m)  H*  (m)  [  Hm)P(m)Hl  (m)  +  N(m)]  ^H(m)P(m)  (34) 

where 

P(m)  =  j^(m,  m-  1)  V(m-  1 )  (m,  m-  1 )  +  G(m)  W(m)  G’(m)  (35) 
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The  most  evident  computational  advantage  in  determining  V(m)  by  Eq.  34 
rather  than  Eq.  31  is  that  Eq.  34  requires  the  inversion  of  a  single  p  Xp 
matrix;  whereas,  Eq.  31  requires  the  inversion  of  three  matrices,  one  of 
order  p  Xp  and  two  of  order  m  X  m.  A  less  evident  advantage  exists  because 
the  number  of  observables  is  very  often  less  than  the  order  of  the  state  vector. 
In  this  instance,  p  <  m  and  Eq.  34  requires  the  inversion  of  a  smaller  matrix 
than  Eq.  31.  The  advantage  of  the  alternative  expression  for  x(m|m),  Eq.  33 
rather  than  Eq.  30,  is  that  the  inverse  of  P(m)  is  not  required. 

Just  as  in  the  previous  examples,  V(m)  does  not  depend  on  the  observed 
data  and  can,  therefore,  be  precomputed.  Then  because  only  simple  matrix 
additions  and  multiplications  are  required,  x(m|m)  can  be  rapidly  computed 
as  new  data  becomes  available.  Furthermore,  the  error  performance  can  be 
investigated  without  carrying  out  a  complete  simulation. 

Collecting  all  results,  we  have 

Model:  x(i)  =  ff(i,  i-  1)  x(i-  1)  +  G(i)  w(i) 

r(i)  =  H(i)x(i)  +  n(i) 

for  i  =  1,  2,  .  .  .  ,  m 

j 

j  Estimation  Equations: 

(1)  x(m  |  m)  «  x(m  |  m  -  1 )  +  V(m)  H’ (m)  N  ^(m)[r(m)  -  H(m)x(m|m-1)] 
where 

x(m  |  m  -  1 )  =  ^(m,  m- 1 )  x(m- 1  |  m- 1 ) 

(2)  V (m)  =  [H'(m)N  1(m)H(m)  +  P  ^(m)]  1 

!  _  .j 

=  P(m)-  P(m)  H'(m)[H(m)  P(m)  Hf(m)+N(m)]  H(m)P(m) 

where 

P(m)  =  ^(m,  m-l)V(m-l)^1  (m,  m  -  1 )  +  G(m)  W(m)  G’ (m) 
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Equations  33  and  34  are  difference  equations  for  the  optimum  estimate  and 

Z  3 

were  first  derived  by  Kalman  and  Bucy  in  1961  *  .  The  associated  initial 

conditions,  x(0|0)  and  V(0),  are  determined  from  the  a  priori  knowledge  of 
x  ,  the  initial  state.  The  known  a  priori  mean  is  x(0|0),  E[x  J  =  x(0):  and 
the  a  priori  covariance  matrix  is  V(0),  Efx^—  xq]  [x^  —  x(0)]  1  =  V(0). 

This  brief  derivation  of  the  discrete  Kalman-Bucy  equations  ignores 
several  important  issues  which  range  from  questions  of  the  existence  of 
inverse  matrices  to  the  uniqueness,  stability,  and  asymptotic  behavior  of 
x(m|m)  and  V(m).  These  issues  are  discussed  in  detail  in  the  cited  references. 
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3.  NONLINEAR  ESTIMATION  THEORY 
3.  1  The  Nonlinear  Estimation  Model 

The  nonlinear  estimation  model  differs  from  the  linear  estimation  model 
only  to  the  extent  of  including  a  possibly  nonlinear  transformation  of  the  state 
vector  in  the  observations.  The  new  model  is  described  by  the  relations: 

x(i)  =  ^(i,  i-  1)  x(i-  1)  +  G(i)  w(i)  (36) 

x(0)  =  xq 

and 

r(i)  =  h [i  :  x(i)]  +  n(i)  (37) 

for  i  =  1,  2,  .  .  .  ,  m.  We  assume  h[i  :  x(i)]  is  suitably  restricted  so  that,  for 
instance,  it  has  a  Taylor  series  expansion 


h(i  :  x)  =  h(i  :  z)  + 


<xk-zk>  sr 


+  higher  order  terms 


x  =  z 


Sh 

=  h(i  :  z)  +  (-^—  )  (x  —  z)  +  higher  order  terms 
ox 
—  z 


(38) 


Sh 

where  (-^—)  is  a  matrix  of  partial  derivatives  of  h(i  :  x)  evaluated  at  the 
ox  —  — 

~  -  Sh 

point  x  =  z;  the  k-row,  jfc-column  element  of  the  matrix,  ,  is  —  h,  (i  :  x). 

—  —  O  X  OX  +  K  — 


This  matrix  is  commonly  referred  to  as  the  Jacobian  matrix. 

3.  2  The  Nonlinear  Estimation  Problem 

Just  as  in  the  linear  case,  we  seek  to  estimate  x(m),  the  state  at  the 
(moving)  endpoint  of  the  observation  interval,  based  on  all  the  accumulated 
observations,  r^  The  procedure  we  shall  take  is  the  following. 
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We  assume  first  that  the  trajectory  traced  out  by  the  sequence  of  points 


x(l),  x(2),  .  .  .  ,  x(m)  does  not  differ  greatly  from  some  trajectory  described, 
for  instance,  by  z(]),  z(2),  .  .  .  ,  z(m).  Then,  for  each  i,  we  expand  h  [i  :  x(i)] 
in  a  Taylor  expansion  about  z(i)  and  retain  only  the  first  two  terms.  The 
observations  are  thereby  linearized  about  the  trajectory  z(l),  z(2),  ...»  z(m), 
and  the  results  of  the  preceding  section  (i.  e.  ,  the  Kalman-Bucy  equations)  can 
be  used  to  estimate  the  state  at  the  end  point  of  the  observation  interval,  x(m). 

The  choice  of  the  trajectory  about  which  the  linearization  is  performed 
depends  on  the  particular  application.  We  can  describe  the  sequence  z(l), 
z(2),  .  .  .  ,  z(m)  as  being  prespecified  or  not  depending  upon  whether  or  not  it 
is  known  before  data  is  received.  The  technique  of  linearizing  about  a  pre¬ 
specified  trajectory  and  then  applying  the  Kalman-Bucy  equations  to  estimate 
the  true  trajectory  has  been  used  with  success  in  a  variety  of  space  applications 
in  which  a  vehicle  was  controlled  so  as  to  follow  a  prescribed  path,  such  as  a 
path  to  the  moon.  *  The  technique  is  also  applicable  to  the  high-altitude  satel¬ 
lite  navigation  system  discussed  below  when  the  user’s  vehicle  is  controlled  so 
as  to  follow  a  prescribed  course.  Such  a  situation  arises,  for  example,  with 
the  point-to-point  travel  of  commercial  aircraft  through  prescribed  air  corridors. 

Very  often,  on  the  other  hand,  no  prespecified  trajectory  is  available 

and  z(i),  for  each  i,  must  be  generated  as  data  is  received.  An  example  of 

this  situation  arises  in  the  navigation  context  when  the  user  does  not  follow  a 

prescribed  course  as  in  a  tactical  or  evasive  maneuver.  If  the  actual  and 

estimated  trajectories  do  not  differ  greatly,  then  it  is  natural  to  consider  a 

*  The  notion  was  introduced  in  this  context  by  Smith  et  al.  [7]  and  McLean 
et  al.  [  8]  . 
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linearization  about  the  estimated  trajectory  itself.  That  is,  for  each  i,  we  can 
let 


z(i)  =  x(i  |  i)  i  =  1,  2,  .  .  .  ,  m  (39a) 

where  x(i  |  i)  is  the  estimate  of  x(i)  based  on  all  the  data  available  at  time  i, 

— 1.  i 

We  shall  see  below  that  this  choice  for  z(i),  Eq.  39a,  leads  to  nonlinear 
equations  for  the  estimate  and  that  these  equations  must  be  solved  simultaneously. 
Since  solving  simultaneous,  nonlinear  equations  is  generally  difficult,  even  with 
the  aid  of  a  computer,  we  shall  propose  an  alternative  choice  for  z(i)  that 
obviates  this  particular  problem.  The  alternative  choice  will  work  nearly  as 
well  as  that  of  (39a)  under  most  circumstances,  certainly  most  circumstances 
where  a  linearization  procedure  will  work  at  all.  Let 


z(i)  =  x(i  |  i-1) 


(39b) 


where  x(i  |  i-1)  is  the  extrapolated  or  predicted  value  of  x(i)  based  on  r.  . 

1 ,  1  -  I 

which  includes  all  of  the  data  available  prior  to  time  i.  The  advantage  of  this 
choice  is  that  while  we  are  still  led  to  nonlinear  equations,  they  can  be  solved 
sequentially  in  a  straightforward  manner. 

The  sequence  z(l),  z(2),  .  .  .  ,  z(m)  can  also  be  viewed  as  the  mean  or 
the  conditional  mean  (i.  e. ,  estimated  mean  of  the  stochastic  sequence  x(  1),  x(2), ...  , 
x(m).  From  this  viewpoint,  a  !,pr  especified"  sequence  z(i)  corresponds  to  the 
case  of  a  known  mean  from  which  it  is  assumed  x(i),  i  =  1,2,...,  m,  does  not 
differ  greatly.  A  non-pre specified  sequence  z(i)  then  corresponds  to  the  case 
of  an  unknown  mean  that  must  be  estimated  as  data  is  accumulated. 
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Let  us  now  illustrate  these  procedures  by  considering  the  nonlinear 
estimation  problem.  First  we  shall  examine  a  special  case  paralleling  the 
linear  estimation  problem  of  Section  Z.  Z.  1;  that  is,  we  shall  remove  all  time 
dependence  and  employ  only  one  observation.  Then  we  shall  reintroduce  time 
dependence  and  treat  the  general  nonlinear  problem. 

3.  Z.  1  Special  Case:  Time-Independent  Observations 
Let 

r  =  h(x)  +  n  (40) 

where  x  and  n  are  Gaussian  random  variables  with  known  means,  x^  and  0, 
and  covariance  matrices,  W  and  N,  respectively.  Our  reason  for  choosing 
this  special  case  is  that  it  models  the  high -altitude  satellite  navigation  system 
discussed  in  Section  IV  when  only  one  observation  is  taken  so  that  changes  in 
the  relative  position  of  the  user  and  satellites  need  not  be  considered.  We 
assume  that  x  docs  not  deviate  significantly  from  its  mean,  x^,  and  then 
expand  h(x)  about  this  point.  To  a  r.l^^p  aDDr^::L nation  we  then  have 

dh 

r  -  h(x  )  =  (■^-)  (x  -  x  )  +  n  (41) 

—  —  — o  o  X  —  — o  — 

—  X 
— o 


We  now  want  to  apply  the  results  of  the  linear  estimation  problem  of 

Section  Z.  Z.  1.  To  this  end,  we  note  that  r  —  h(x  )  may  be  taken  as  the 

dh  -  -  -o 

observed  signal,  as  the  matrix  H  =  H(x  ),  and  x  —  x  as  a  zero  mean 

&  dx  x  —  —  — o'  -  — o 

- o 

Gaussian  random  variable.  It  then  follows  from  Eq.  13  that 

dh 

x*  =  V*(x  )  (— )  N_i  [r  -  h(x  )]  (4Z) 

— *  —  — o  OX  —  — - o 

—  X 

— o 


19 


where  the  asterisk  on  x*  indicates  that  the  estimate  is  an  approximation  to 
the  optimum  estimate,  x;  we  shall  refer  to  x*  as  being  "quasi -optimum.  " 
The  matrix,  V*(x^),  is  given  from  Eq.  1Z  by 


_  1  dh  _  ,  dh  1 

i*<5o»  *  i r  +  y  r  <3?  r‘ 

—  X  —  X 

— o  — o 


(43) 


An  alternative  expression  for  V*(x^)  can  be  obtained  by  using  Eq.  14. 

The  error  covariance  matrix  V,  defined  by 

V  =  E [x  -  x*]  [x  -  x*]  '  (44) 

equals,  or  nearly  equals,  V*(x^)  when  the  error,  x—  x*,  is  small.  ^  More 
generally,  the  Cramer-Rao  bound  ^  can  be  used  to  demonstrate  that  V  satisfies 
the  inequality 

V  £  V*(x  )  (45) 

—  —  — o 

whatever  the  size  of  the  error.  This  inequality  implies  that  the  error  ellip¬ 
soids  associated  with  V, 

[x  —  x*]  1  V  [x  —  x*]  =  constant,  (46) 

never  lie  within  those  associated  with  V*(x  ) 

—  — o 

[x  —  x*]  1  V*  (x^)[x  —  x*]  =  constant  (47) 

for  the  same  constant.  It  follows  that  V*(x^)  can  be  used  to  study  the  best 

attainable  performance  of  x*  as  an  estimate  of  x.  For  this  purpose,  various 

norms  of  V *(x  )  can  be  used.  Two  possible  norms,  which  were  mentioned  in 
—  — o 

+  We  say  x  —  x*  is  small  if  some  related  norm  is  small.  A  typical  norm 
which  can  be  used  is  tr[x  —  x*]  [x  —  x*]  *. 


ZO 


the  linear  case  of  Section  2.  2.  1,  are  the  largest  eigenvalue  and  trace  of 

— *  (— O  )  ' 

3.  2.  2  General  Case:  Time  -  Dependent  Observations 

For  this  case,  we  have  the  observed  sequence  given  by  Eq.  37  as 

r(i)  =  h[i :  x(i)]  +  n(i)  (37) 

for  i  =  1,2,...,  m.  We  now  expand  h  [  i  :  x(i)]  about  z(i),  for  each  i,  and 
then  use  the  small  error  assumption.  At  this  point,  z(i)  can  be  either  pre¬ 
specified  or  not.  We  shall  consider  the  separate  cases  below.  To  a  close 
approximation,  the  result  is 

.  3h  3h 

P(i)  =  r(i)  -  h[i:z(i)]  +  (5=)  z(i)  =  (^=)  x(i)  +  n(i)  (48) 

-  z(i)  -  z(i) 


for  i  =  1,2,...,  m.  The  results  of  the  general  linear  estimation  problem  can 
be  applied  when  we  identify  p(i)  as  the  observed  signal  and  (&h/&x)  ...  as  the 

—  —  —  z^ 1 ) 

matrix,  H(i).  From  Eq.  33  we  then  obtain 


3h  ,  ,  3h 

x*(m|m)  -  x*(m|m-l)  +  V*(m)  (— )  N  (m)  {p(m)  —  (— )  x*(m|m-l)} 

“  z(m)~  ”  a-z(m)“ 


(49) 


=  x*(m  |  m  -1) 

3h  ,  i  3h 

+  V*(m)(— )  N~  (m){r(m)  —  h[m:z(m)]  —  (-^— ) 

z(m)“  -  -  -  Sxz(m) 


X  [x*(m  I  m-1)  —  z(m)]  } 


(50) 


where  V*(m)  is  a  matrix  specified  from  Eq.  3 1  by 


V*(m)  = 


1  i  3h 

—  z(m)  —  z(m) 


-  1 


(51) 
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and  where 


P(m)  =  ^(m,  m-  1 )  V*(m-  1 )  (m,  m -  1 )  +  G(m)  W(m)  G'(m). 

An  alternative  expression  for  V*(m)  can  be  given  by  using  Eq.  34.  A  further 
simplification  of  Eq.  50  is  possible  if  it  is  assumed  that  z(m)  does  not  differ 
greatly  from  the  extrapolated  state,  x*(m  |  m-1): 

z(m)  ^  x*(m|  m-1)  ( 5 Z ) 


Then 

Sh 

h[m:  x*(m|m-l)]  «  h[m:z(m)]  +  )  f  x*(m  |  m  -1)  —  z(m)] 

-  z(m)  " 

and  the  expression  for  x*(m  |  m)  becomes 


Sh  , 

|  m)  =  x*(m  |  m-1)  +  V*(m)  N  (m)  {  r(m)  —  h  f  m: x*(m  |  m -1)]  } 

-  z(m)  ~ 

Equations  51  through  53,  along  with  the  initial  conditions, 

x*(0|0)  =  E[x(0)]  =  x  (54) 


and 

V*(0)  =  E[x(0)  —  xQ]  [x(0)  —  Xq]  1  4  V(0)  (55) 


specify  the  qua  si -optimum  estimate  of  the  state  at  the  endpoint  of  the  observa¬ 
tion  interval.  It  can  be  observed  now  that  the  trajectory  followed  by  z(i),  for 

each  i,  enters  into  the  estimate  through  the  Jacobian  matrix,  (gh/gx) 

-  —  z(m) 

which  occurs  both  in  Eqs.  51  and  53. 

In  the  case  of  prespecified  trajectories,  z(i),  for  each  i,  is  known  in 

advance  so  that  (gh/gx)  ,  *  and,  therefore,  V*(m)  can  be  precomputed.  We 

—  —  /  — 

then  observe  that 


(53) 
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(1)  x*(m  |  m)  can  be  determined  rapidly  and  simply  as  new  data 

becomes  available.  This  is  because  the  "gain"  matrix,  V*(m)(dh/dx)  N"l(m), 

z(m) 

can  be  precalculated  so  that  the  only  "on-line"  matrix  operations  required  in 
Eq.  53  are  additions  and  multiplications  and  not  inversions. 

(2)  The  error  performance  can  be  investigated  using  V*(m)  without 
making  any  observations  or  involved  computer  simulations. 

On  the  other  hand,  V*(m)  cannot  be  precomputed  when  there  is  no  pre¬ 
specified  trajectory;  that  is,  when  z(i)  must  be  generated  as  data  is  received. 

If  z(i)  is  defined  by  (39a): 

z(i)  =  x*(i  |  i)  i  =  1,  2,  .  .  .  ,  m  (Eq.  39a,  repeated) 

then  (51)  and  (53)  are  coupled  and  must  be  solved  simultaneously  Alterna¬ 
tively,  if  z(i)  is  defined  by  (39b), 

z(i)  =  x*(i  I  i-1)  i  =  1,  2,  .  .  .  ,  m  (Eq  39b,  repeated) 

then  (51)  and  (53)  are  coupled  but  they  can  be  solved  sequentially.  The  error 
performance  cannot  be  investigated  in  either  of  these  cases  without  an  involved 
computer  simulation  Note  also  that  if  x*(m  |  m)  does  differ  from  x*(m  m-1) 
enough  so  that  (9h/3x)^^  \  ^  is  appreciably  different  from 

then  solving  Eqs.  (51)  and  (53)  sequentially  is  not  sufficient.  The  linearization 
of  h  about  the  final  estimate  must  agree  with  the  linearization  used  in  the  solu¬ 
tion  of  these  equations. 
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4.  APPLICATION  OF  NONLINEAR  ESTIMATION  THEORY  TO  THE  PROBI  EM 


OF  NAVIGATING  WITH  HIGH- ALTITUDE  SATELLITES 

We  now  want  to  apply  the  results  derived  in  the  preceding  sections  to 
simplified  versions  of  the  high-altitude  satellite  navigation  system  described 
in  the  introduction.  We  shall  study  the  problem  in  two  setps.  First,  we 
examine  the  simplest  case  where  only  one  observation  is  taken  so  that  user 
and  satellite  position  changes  need  not  be  considered.  Then,  second,  we  shall 
introduce  multiple  observations  and  include  simple  user  motion. 

4.  1  Special  Case  I:  Fixed- Position  and  Single  Observation 

The  geometry  associated  with  the  navigation  system  for  the  case  of 
three  fixed  satellites,  one  fixed  control  station,  and  one  fixed  user  is  shown 
in  Fig.  1.  We  shall  use  an  earth-centered  rectangular  coordinate  system. 

The  position  vectors  associated  with  the  satellites,  control  station,  and  user 
are  defined  in  Table  I. 


Position 

Satellites 

j=  1,  2,3 

Control  Station 

He 

User 

<  _ 

Hu 

Table  I 


Let  h  be  the  user's  height  measured  from  the  earth's  center: 

1/2 


h  =  iRu 


=  lE'uEj 


(56) 


h  is  nominally  equal  to  the  radius  of  the  earth. 
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Figure  1. 


NAVIGATION  SATELLITES 


Geometry  of  the  high-altitude  satellite  navigation  system 


-25- 


The  user-to-satellite  and  control  -  station-to- satellite  distances  indi¬ 
cated  in  Fig.  1  are  given  by 

sj  =  1  p j  —  p u  I  *  f(Ej  "Ej'tej  “Eu^1^2  (5?) 

and 

lj  =  lEj-£j  =  ^Ej-Ec),(Ej-Eu)}l/2  (58> 

for  j  =  1,  Z,  3,  respectively. 

th 

Let  t.  be  the  time  difference  between  the  clock  of  the  j  satellite  and 


a  master  clock.  The  vector 


T  = 


(59) 


is  then  a  time-error  vector  associated  with  the  three  satellite  clocks.  Simi¬ 
larly,  let  tq  be  the  time  error  associated  with  the  user's  clock  (again  rela¬ 
tive  to  the  master  clock).  We  assume  that  for  i  =  0,  1,  2,  and  3,  is 
measured  in  distance  units  [see  footnote  (*)  on  next  page]  . 

The  unknown  parameter  (or  state)  vector  describing  the  navigation 
system  is  given  in  terms  of  the  defined  quantities  as: 


x  = 


(19  components) 


(60) 
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The  user  attempts  to  determine  his  position  by  observing  the  control 
station  monitoring  signals,  the  satellite  timing  signals  (relative  to  his  clock), 
and  his  height.  We  shall  simplify  the  problem  at  this  point  by  assuming  that 
the  control  station  can  make  very  accurate--in  fact,  perfect- -measurements 
of  the  satellite  positions  and  satellite  clock  errors.  In  this  instance,  the 
control  station  positions,  p^,  and  satellite  clock  errors,  t»  are  eliminated 
from  the  problem  and  the  user’s  observable  may  be  taken  as 

r  =  h(x)  +  n  (61) 


where 


— 

— 

J— 

- 

h^x) 

S1 

+ 

T 

O 

1  Hi 

-Eu 

|  + 

T 

O 

h2(x) 

S2 

+ 

T 

O 

k2 

-Eul 

+ 

T 

O 

h3(x) 

S3 

+ 

To 

Ie3 

-Eu' 

'  + 

T 

O 

h4(x) 

h 

Ieu' 

and  n  is  the  observation  noise  with  an  assumed  covariance 


E[nn'] 


0 


Here  a  is  the  ranging-error  variance  that  results,  for  example,  from 
receiver  noise,  propagation  effects,  and  short-term  clock  instabilities; 
is  the  height  error  variance. 

The  vector  x,  defined  by: 


*  (from  previous  page)  Just  as  radar  echo  delays  are  converted  to  distances 
by  using  the  speed  of  light. 
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X  = 


Pu 


Hi 


Hz 

£3 


(13  components) 


is  assumed  to  have  a  mean  E(x)  =  Xq  which  accounts  for  the  nominal  posi¬ 
tions  of  the  satellite  and  user  and  also  the  nominal  value  of  the  user's  clock 
error.  The  a  priori  error  covariance  matrix  is  defined  by 

W  =  E(x  —  x^)(x  —  Xq)  1  (dimension  13  X  13) 

The  diagonal  components  of  W  are  the  a  priori  variances  associated  with  the 
satellite  and  user  positions  and  with  the  user's  clock. 

Define  the  Jacobian  matrix  by 


bhl 

dhj 

8h 

8x2 

&x 

8h 

Sh2 

Sh2 

8h 

dx 

dx^ 

8x2 

dx 

8h4 

8h4 

8h 

8x^ 

8x2 

bx 

13 


13 


13 


(62) 


Then,  from  (43),  V*(x  ),  the  error  covariance  matrix  associated  with  esti¬ 


mating  x  when  E[x]  =  x  ,  is  given  by 
—  —  — o 


V*(x0)  = 


/  _CL 


-1 


X 


lbh\ 

/Shi 

/ 

-1  (bh\ 

I  _ 

w 

+  N 

_ 

1  8x 

„  — 

dx 

\  -1 

*0  '  "I 

IX 

o 

\  *1 

=  w  -  w 


w 


dx 


-0 


(63) 
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Equation  (63)  is  being  used  in  a  computer  study  of  the  navigation  problem 


for  the  special  case  described.  Various  user  positions,  x^,  and  a  priori 
statistics,  W,  are  being  examined.  The  results  will  appear  in  a  companion 
report. 

4.  Z  Special  Case  II:  Simple  Motion  and  Multiple  Observations 

In  this  example,  we  shall  make  the  following  assumptions: 

(i)  The  control  station  can  make  very  accurate  measurements  of 
the  satellite  positions  and  clock  errors;  this  removes  p^  and  t  from  the 
problem. 


(ii)  There  are  three  fixed  satellites. 

(iii)  The  user  undergoes  motion  along  some  nominal  prescribed 
trajectory;  that  is,  his  mean  trajectory  is  known.  His  "state,  M 


s 

— u 


(i) 


(i) 


v 

— u 


i  =  1 ,  Z,  .  .  .  ,  m 


(64) 


consisting  of  his  position  and  velocity,  satisfies 

s  (i)  =  $  (i,  i-1)  s  (i-1)  +  G  (i)  w  (i)  i  =  1,  2,  .  .  .  ,  m  (65) 

— u  — u  — u  — u  — u 

where  $^(i,  i - 1 ) *  G^(i),  and  w^(i)  are  the  sarne  as  defined  by  Eq.  (1)  of 
Section  Z.  1  except  that  here  we  allow  w^(i)  to  have  a  nonzero  mean.  The  non¬ 
zero  mean  of  w  (i)  is  such  that  s^(i),  i  =  1,  2,  .  .  .  ,  m,  follows  the  prescribed 

trajectory  on  the  average;  that  is,  we  are  assuming  that  the  (known)  mean  of 

the  stochastic  sequence  s^(i),  i  =  1,  Z . m,  is  the  prescribed  trajectory  about 

which  the  linearization  is  carried  out. 

(iv)  The  user's  clock  is  stable  but  has  a  fixed  offset,  t  ,  from 
'  o 

true  time  (as  indicated  on  a  master  clock). 
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(v)  Multiple  observations  are  taken. 

The  state  vector  x(i),  i  =  1,  Z,  .  .  .  ,  m,  associated  with  the  estimation 
problem  is  defined  by 

P 1  (i) 


x(i)  = 


P2(i) 

P3(i) 


(16  components) 


(66) 


where,  because  of  the  above  assumptions, 


I 

0 

0 

0 

0 

0 

0 

I 

0 

0 

0 

0 

x(i)  = 

0 

0 

I 

0 

0 

x(i-l)  + 

0 

0 

0 

0 

gu(i.i-l) 

0 

G  (i) 

— u 

0 

0 

0 

0 

1 

0 

_ 

_ 

w  (i) 
— u 


=  0(i,  i-l)x(i-l)  +  G(i)  w^li)  for  i  =  1,  2 . m  (67) 


The  observed  sequence  satisfies 
r(i)  =  h[x(i)]  +  n(i) 


i  =  1,2 . m 


(68) 


where 


h  f  x(i)  ]  = 


El  -  Eu^  I  +  T0 

p2  -  pu(i)  I  +  Tq 

p3  -  Eu(i>  I  +  To 
I  Pu(i)  1 


(69) 
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and  n(i)  is  the  observation  noise  with  an  assumed  covariance 

2 


E[n(i)n'(j)]  = 


2  0 


0  ar 


6. . 

ij 


(70) 


We  have  assumed  for  this  simple  model  that  the  additive  observation  noise  is 
an  uncorrelated  sequence.  In  practice,  correlated  sequences  would  also  be 
encountered  and  these  would  be  treated  by  expanding  the  dimension  of  x(i). 
Define  the  Jacobian  matrix  by 


ah 

ahj 

Bh 

hxi 

Sxz 

Sx 

*h2 

Sh2 

Sh 

9x  j 

$x2 

bx 

sh4 

ah 

Sx2 

dx 

13 

2 


13 


13 


(71) 


x(i) 


Then,  from  (51),  the  error  covariance  matrix  associated  with  estimating  x(m) 
when  E[x(m)]  =  z(m),  is  given  by 


V*(m)  = 


P"  1  (m)  + 


/a  h\' 


N"  !(m) 


/ah' 

1  dx 

z(m) 

\  -( 

-1 


(72) 


where 


P(m)  —  $(m,  m-1)  V  *  ( m  -  1 )  §  1  ( m ,  m-1)  4-  G(m)  W(m)  G’(m) 

and  W(m)  is  the  a  priori  covariance  matrix  defined  by 
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W(m)  =  E  [  w(m)  —  E  w ( m ) ]  [  w ( m )  —  E  w ( m ) ] *  1 2 3 

Note  that  V*(m)  is  determined  from  the  difference  equation  (72)  by  starting 
with  the  initial  condition  V*(0)  =  Vq,  the  a  priori  error  covariance  matrix,  and 
serially  calculating  V*(l),  V*(2),  .  .  .  ,  V*(m). 

Equation  (72)  is  being  used  in  a  computer  study  of  the  navigation  prob¬ 
lem  for  the  special  case  described.  Various  simple  user  trajectories  and 
a  priori  statistics,  w(m),  are  being  examined.  The  results  will  appear  in  a 
companion  report. 

Preliminary  Conclusions 

We  have  indicated  with  two  simple  examples  how  the  nonlinear  estima¬ 
tion  procedure  of  section  3  can  be  used  for  the  problem  of  navigating  with 
high-altitude  satellites.  Several  questions  that  must  be  addressed  before  the 
technique  can  be  applied  successfully  to  the  more  general  navigation  problem 
are: 

(1)  Models  must  be  developed  for  describing  the  motion  of  satel¬ 
lites  placed  in  a  near  -  synchronous  orbit  and  for  characterizing  the  errors 
made  in  determining  their  coordinates; 

(2)  Models  must  be  developed  for  statistically  describing  the 
effects  of  timing  errors  due  to  atmospheric  effects,  such  as  refraction.  These 
effects  generally  depend  on  the  relative  user -to -  satellite  position.  Conse¬ 
quently,  the  associated  covariance  matrices  will  depend  on  the  state  (x); 

(3)  An  investigation  of  the  sensitivity  of  the  navigation  model  to 
changes  in  assumptions  or  parameters  (such  as  the  a  priori  error  covariance 
matrix)  should  be  made. 
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Appendix 


We  wish  to  establish  the  lemma: 


If 


where 


s'1  =  T"1  +  M' u'1  M  (Al) 

S  \  T  \  and  U  *  are  symmetric  po sitive -definite  matrices;  then 
S  =  T  -  TM'[M  TM'  +  U]  '  1  M  T  (AZ) 


Proof: 

The  existence  of  the  inverse  matrices  required  for  Eq.  A 2  follows  from 
the  po  sitive -definitene  s  s  of  S  T  and  U  ^ .  We  establish  the  Lemma  by 
direct  substitution. 

S  S'  1  =  [  T  -  T  M'(M  T  M'  +  U)' 1  M  T]  [  T' 1  +  M' u' 1  M] 

=  J_-  TM'[(MTM-  +  U)'1  M  -  U"1  M  +  (MTM1  +  U)"  1  M  T  M'  U'1  M] 
=  TM'HMTM1  +  U)'1  U  -  J_+  (MTM1  [(U)'lMTM']  U"1  M 
=  J_ -  TM'[(MTM'  +  U)'1  (MTM1  +  U)  -  1]  U'  1  M 
=  I 
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